Coarse-grained analysis of a lattice Boltzmann model for planar streamer fronts 
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We study the traveling wave solutions of a lattice Boltzmann model for the planar streamer fronts 
that appear in the transport of electrons through a gas in a strong electrical field. To mimic the 
physical properties of the impact ionization reaction, we introduce a reaction matrix containing 
reaction rates that depend on the electron velocities. Via a Chapman-Enskog expansion, one is able 
to find only a rough approximation for a macroscopic evolution law that describes the traveling wave 
solution. We propose to compute these solutions with the help of a coarse-grained time-stepper, 
which is an effective evolution law for the macroscopic fields that only uses appropriately initialized 
simulations of the lattice Boltzmann model over short time intervals. The traveling wave solution is 
found as a fixed point of the sequential application of the coarse-grained time-stepper and a shift- 
back operator. The fixed point is then computed with a Newton-Krylov Solver. We compare the 
resulting solutions with those of the approximate PDE model, and propose a method to find the 
minimal physical wave speed. 



I. INTRODUCTION 

When a gas of neutral atoms or molecules is exposed 
to a strong electrical field, a small initial seed of electrons 
can lead to an ionization avalanche. Indeed, the seed elec- 
trons are accelerated by the field and gain enough energy 
to ionize the neutral atoms when they collide. The two 
slow electrons that emerge from this reaction, i.e. the 
impact and the ionized electron, are again accelerated by 
the field and cause, on their turn, an ionization reaction. 
Simultaneously the electrical field is locally modified be- 
cause of the charge creation. This interplay between the 
dynamics of the electrons and the electrical field can lead 
to a multitude of phenomena studied in plasma physics 
such as arcs, glows, sparks and streamers. 

In this article we will focus on the initial field driven 
ionization that can lead to traveling waves known as 
streamer fronts. These waves have previously been stud- 
ied by Ebert et al. [l[ who introduced and analyzed the 
minimal streamer model, a one-dimensional model for 
the propagation of planar streamer fronts. This model 
consists of two coupled non- linear PDEs: a reaction- 
convection-diffusion equation for the evolution of the 
electron density and a Poisson-like evolution equation 
for the electrical field. The reaction term is based on 
the Townsend approximation that expresses the growth 
of the number of electrons as a function of the local elec- 
trical field. 

During the last two decades, however, a lot of progress 
has been made in the microscopic understanding of im- 
pact ionization reactions in atomic and molecular sys- 
tems. In this reaction an impact electron ionizes the tar- 
get and kicks out an additional electron. There are sev- 
eral successful theories that can predict the exactprob- 
ability distribution of the escaping electron 0, H, [§, Q3 • 
In the next decade, we expect that the theoretical tools 
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will be able to accurately predict the microscopic physics 
of electron impact on molecular targets such as N2 and 
O2, the most important molecules in the composition of 
air. This progress in the understanding of the impact 
ionization reaction, however, has not been incorporated 
in the description of the macroscopic behavior such as the 
minimal streamer front of Ebert et al. Instead, such mod- 
els still make use of a phenomenological approximation 
to the reactions, such as the Townsend approximation. 
This article extends the minimal streamer model and in- 
corporates more microscopic information. We model the 
system by a Boltzmann equation, which is constructed 
such that the cross sections in the collision integral re- 
semble the true microscopic cross sections. 

To find the traveling wave solutions of this more micro- 
scopic model, we exploit a separation of time scales be- 
tween the relaxation of the electron distribution function 
to a local equilibrium and the evolution of the macro- 
scopic fields (electron density and electrical field). It is 
known from kinetic theory that the first process is fast: 
once initialized, it takes a molecular gas not more than 
a few collisions to relax to its equilibrium state. 

In kinetic theory, the fast times scales are often elimi- 
nated from the problem by assuming a local equilibrium 
distribution function which leads to a reaction-diffusion 
model with transport coefficients that depend on the lo- 
cal electrical field. This reduction method, however, is 
only successful in the absence of steep gradients in the 
electron density [TH, an assumption not valid for the 
planar streamer fronts that have, typically, very steep 
increases in the electron density. 

In this article, we take an alternative route and find 
the traveling wave solution through a so-called coarse- 
grained time-stepper (CGTS) that exploits the separation 
of time scales to extract the effective macroscopic behav- 
ior. This method was proposed by Kevrekidis et al. || 
and the numerical aspects of its application to find trav- 
eling wave solutions of lattice Boltzmann models have 
recently been studied ||. The time-stepper uses a se- 
quence of computational steps to evolve the macroscopic 
state. This sequence involves: (1) a lifting step, which 
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creates an appropriate electron distribution function for 
a given electron density, (2) a simulation step, where 
the lattice Boltzmann model is evolved over a coarse- 
grained time step AT, and (3) a restriction step, where 
the macroscopic state is extracted from the electron dis- 
tribution function. This method does not derive effec- 
tive equations explicitly, and therefore allows steep gra- 
dients to be present. We compare our results with those 
obtained by deriving an approximate macroscopic PDE 
model through the more traditional Chapman-Enskog ex- 
pansion. The paper therefore illustrates the applicability 
of the coarse-grained time-stepper on a non-trivial prob- 
lem where the exact macroscopic equations are hard to 
derive. 

The outline of the paper is as follows. In section II, we 
shortly review the physics of the impact ionization reac- 
tion and the streamer fronts, and recapitulate the Boltz- 
mann equations. In section III, we derive from the Boltz- 
mann equation a lattice Boltzmann model with multiple 
velocities and discuss how the ionization reaction, exter- 
nal field and electron diffusion are incorporated in the 
model. Section IV derives a macroscopic PDE from the 
model using the Chapman-Enskog expansion and dis- 
cusses the minimal velocity of the traveling waves. Sec- 
tion V formulates the coarse-grained time-stepper and 
VI how the traveling wave solutions arc found. Finally 
in section VII, we have some numerical results. 



II. MODEL 
A. The physics of the impact ionization reaction 

The impact ionization reaction is a microscopic reac- 
tion where electrons with, typically, an energy around 
50eV collide with an atom or a molecule and ionize this 
target. The reaction rates of this process depend sensi- 
tively on a number of parameters. Let us consider the 
simplest system: an electron hits a hydrogen atom with 
a bound electron in its ground state. When the incoming 
electron has an energy larger than the binding energy of 
the electron in the atom, it can kick out, with a certain 
probability, the bound electron that will escape, together 
with the impact electron, from the atom. The total en- 
ergy of the two electrons after the collision is equal to 
the energy of the incoming electron minus the original 
binding energy of the bound electron. 

The reaction rates of this process are expressed by cross 
sections; these are probabilities that a certain event will 
take place. One such cross section is the triple differ- 
ential cross section, which is the probability to find af- 
ter the collision one electron escaping in the direction 
(6*i, 0i ) with an energy E\ and a second electron escap- 
ing in the direction (Bi,<\>-i) with an energy E 2 , where 
the angles are measured with respect to the axis defined 
by the momentum of the incoming electron. Since the 
two electrons repel each other ,it is more likely that they 
escape in opposite directions (l3j |. 
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Figure 1: A sketch of the typical shape of the impact ion- 
ization cross section ( see the experimental results in [(|). 
On top, we show the total cross section as a function of the 
impact energy where below a threshold energy of 25eV no 
reactions take place. In the proposed model we distinguish 
between slow particles with an energy below this threshold 
that do not react and fast particles with their energy above 
this threshold. The fast reacting particles experience a cross 
section of R, as indicated by the dashed line. In the bottom 
figure, we show the energy differential cross section for the 
escaping electron, where the total energy of the escaping elec- 
trons is 25eV. Since the two electrons are indistinguishable 
there is a symmetry. In our five speed model, we make the 
approximation that the two electrons can only escape with 
equal energy sharing. 



When this cross section is integrated over all angles 
(0i, 4>\) and (#2, 02 ) of the escaping electrons and all pos- 
sible ratios of E\ / E 2 of the electron energies, we get the 
total cross section. This is the total probability that the 
incoming electron will cause an ionization event. This 
total cross section depends on the energy of the incom- 
ing electron and is zero when the energy of the incoming 
electron is below the binding energy of the bound elec- 
tron. Just above this binding energy, there is a steep 
rise in the cross section that is known as a threshold. 
Just above this threshold the cross section is the largest 
and as we further increase the energy the cross section 
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diminishes. This is illustrated in figure [I] (top) . 

When the cross section is integrated of the angles only, 
but not over the relative energies, we get the so called en- 
ergy differential cross section, which is the probability of 
causing an ionization event with a given relative energy 
of the two electrons. In contrast with the electron direc- 
tions, there is no pronounced preference for the energy 
sharing between the two electrons, see figure [T] (bottom) . 
It is only slightly more likely that two electrons will come 
out with unequal energy. 

Recently, several theoretical methods have successfully 
predicted the directions of the escaping electrons, the to- 
tal cross section and the energy differential cross section 
in the hydrogen atom. Wc name exterior complex scal- 
ing , time dependent close cow iling d, HRW-SOW @ 
and convergent close coupling (lOj. 

When the electron hits a molecular system instead of 
an atom, the physics is complicated by the extra degrees 
of freedom. The cross sections now depend on both the 
orientation and internuclear coordinates of the molecule 
at the moment of electron impact, as is seen in processes 
where two electrons are ejected from molecules after it is 
hit by a photon [ll|, 0|. Therefore, there will be some 
random terms in the reaction cross section, which will 
need to be included in a realistic microscopic model. In 
this paper, we will model the microscopic interactions 
using a Boltzmann equation, which is still deterministic; 
extensions that accurately account for random effects will 
be treated in future work. However, we note that the 
coarse-grained time-stepper approach that is used in this 
work has already been applied successfully to study sys- 
tems with stochastic effects 12011 . 



B. Review of the physics of streamer fronts. 

Ebert et al. [l[ introduced the minimal streamer model. 
It consists of two coupled non-linear PDEs: a reaction- 
convection-diffusion equation for the evolution of the 
electron density and an equation that relates the change 
in the electrical field to the charge flux. The electron den- 
sity evolves because of the drift due to the electrical field, 
the electron diffusion and the ionization reaction, which 
is formulated in the Townsend approximation. The re- 
action rate is then given by an exponential that depends 
on the strength of the local electrical field. The evolu- 
tion of the electrical field is determined by Poisson's law 
of electrostatics where the field changes because of the 
charge creation by the ionization reaction. This minimal 
streamer model exhibits both negatively and positively 
charged fronts. The first moves in the direction of the 
electrical field, while the positively charged moves in the 
opposite direction and can only propagate because of the 
electron diffusion and the ionization reactions. Each of 
these fronts appears as a one-parameter family of uni- 
formly translating solutions (since any translate of the 
wave is also a solution). 

In our extension of the Ebert model, we replace the 



reaction-diffusion equation for the evolution of the elec- 
tron density with a Boltzmann equation for the one- 
electron distribution function /(x, v,i), that counts the 
number of electrons in the phase-space volume element 
bounded by position x and x + dx and by speed v and 
v + dv. The Boltzmann equation is 



■v M J y +E(x,t)- ' ' = fi(x,t), 



dv 



0/(x,v,t) 
dt 

(1) 

where E(x, t) is the external electrical field and fi(x, i) 
is the collision operator, an integral operator that inte- 
grates the cross sections of the ionization reaction over 
the velocity space. This Boltzmann equation is coupled 
to an evolution equation for the electrical field. Because 
additional electrons are created, the local charge density 
changes the electrical field through the Poisson law 

V-E(x,t) = 9 (x,t), 

where q(x, t) = (n+ — n e )e/qo is the charge distribution. 
Here, n + represents the number of ions, n e is the number 
of electrons and qo a unit of charge. 

We now connect the change in electrical field with the 
change in the charge distribution. We have 



dg(x,t) 
dt 



V-j(x,i) = 



in which j(x, t) is the charge flux. This leads to an equa- 
tion for the evolution of the electrical field, 



9E(x,i) 
dt 



+ j(x,t) = 0. 



(2) 



Since we assume that the ions are immobile, the flux 
j(x, t) is solely determined by the one-particle distribu- 
tion function /(x, v, t) of the electrons. 

Our extension of the minimal streamer model is now 
the coupled evolution of eq. (JTJ) and (0). Note that the 
set of coupled equations is very similar to the Wigner- 
Poisson problem [THj used to model electron transport 
through diodes. 



III. LATTICE BOLTZMANN DISCRETIZATION 



Together with the impact ionization cross sections, the 
coupled equations ([T]) and @ are a non-linear integro- 
differential equation coupled to a scalar partial differen- 
tial equation for the electrical field. This equation in its 
full dimension is hard to solve, both analytically and nu- 
merically. As a first step, we look at the one-dimensional 
streamer fronts of the Boltzmann equation in the lattice 
Boltzmann discretization. 
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A. Discretization of the Boltzmann equation 

In this section, we discretize the one-dimensional 
Boltzmann equation ([1]). The distribution functions 
f(x,v,t) are discretized on a lattice in space, velocity 
and time. The grid spacing is Aa; in space and At in 
time. The velocity grid of i>i is chosen such that the dis- 
tance traveled in a single time step, ViAt, is a multiple 
of the grid distance Ax, or in short 



.Ax 

'A? 



with 



i G S. 



Typically, only a small set of discrete velocities is used. 
A discretization with three grid points on the velocity 
grid has a set S = { — 1,0,1} and is called a D1Q3 
model. A discretization with five grid points has S = 
{—2, —1, 0, 1, 2} and is a D1Q5 model. The size of the set 
S is denoted by m. We will also denote Cj = ViAt/ Ax 
for the dimcnsionless velocity. 

Note that, for ease of notation, we will also use the set 
S to index matrices. For example, the result of a linear 
operator A working on a vector Vj^s will be denoted as 
SieS AijVj) where both the indices i and j are in S. 
The means that the A-2,-2 matrix element with S = 
{—2, —1,0, 1,2} is the matrix element in the upper left 
corner of the matrix. 



We start from the continuous equation for the distribu- 
tion function in the discrete point (x + At, v{) in phase 
space at time t. In the absence of external forces, the 
Boltzmann equation in this point reads 



df(x + ViAt,Vi,t) df(x + ViAt,Vi,t) 

\-Vi — 



at 

= Q(x + ViAt, Vi, t) 



Ox 



(3) 



A discrete lattice Boltzmann equation is now obtained 
by replacing the time derivative with an explicit forward 
difference, the introduction of an upwind discretization 
of the convection term and a downwind discretization of 
the collision term f2(x + ViAt,Vi,t) and replace it with 
Q(x,Vi,t) 0, 

f(x + Vi At, vi, t + At)- f(x + Vi At, Vi , t) 
At 

+ ^ f(x + «, At, Vj ,t) - f(x, Vj, t) = 
ViAt 

Note that this discretization of the spatial derivative 
becomes less accurate for the largest speeds in the set 
S. Indeed, in the five speed model, for example, the 
largest speed is v± 2 = ±2 Ax /At, and the convec- 
tion term will be calculated from the difference between 
f(x + v± 2 At, v±2,t) and f(x, v± 2 , t), which is 2Ax apart. 
The discretization error is then proportional to 2Ax. 

Equation (|4]) reduces to 



where we have introduced the shorthand fi(x,t) for 
f(x,Vi,t) and Qi(x,t) for Q(x,Vi,t) with i G S. 



B. The collision term 

The collision term consists out of two parts 

Atfli = nf lS + ^ caction (6) 

where the first term will model the electron diffusion, the 
second term the ionization reactions and the third term 
the influence of the external force. Note that we also 
incorporated At in the notation. We will now discuss 
the two terms individually. 

The first term models the electron diffusion as a BGK 
relaxation process [2~H |. In this approximation, it as- 
sumed that the distribution is attracted to a local equi- 
librium distribution function 



1 



nr = — (fi 

T 



(7) 



with f^ q (x,t) the equilibrium distribution for electron 
diffusion. In the five speed model, we choose 



f 

•I t 



eq 

w i p 



with 9 = {0, 1/4, 2/4, 1/4, 0} (8) 



and p(x,t) the electron density. This choice of equi- 
librium weights conserves the number of electrons, but 
does not conserve momentum as a traditional fluid would 
do. Indeed, the electrons diffuse because they randomly 
change their direction during the elastic collisions with 
the much heavier neutral molecular particles. We fur- 
ther chose the weights such that there are no fast parti- 
cles under diffusive equilibrium. The relaxation time r is 
related to the electron diffusion coefficient 



D 



At 



E 



(9) 



Note that, in the literature, the relaxation time r is often 
characterized by its inverse u = 1/r. 

The second term in © is the reaction term J7j eactlon 
that is modelled with a m x m matrix 1Z 



n? a * tion = AtJ2K ij f j , 



(10) 



which represents the velocity dependent reaction rates 
and allows us to select between slow and fast particles. 
For the five speed model, we choose a reaction matrix 



n 



( -R 

R R 



R R 

V -rJ 



(ii) 



fi(x + ViAt, t + At) — fi(x, t) = At Vti(x, t), (5) that describes how the reaction cross sections depend on 
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the velocities of the particles. Each time step At, a frac- 
tion R of the particles with speed v±2 will collide and 
cause a ionization reaction. The reaction rate R is cho- 
sen to match the height of the cross section, see figure 
[TJ Since the colliding fast particles transfer their energy 
to the bound electron, they will loose energy. There- 
fore, the number of particles with speed t>± 2 diminishes 
with a rate — RAt and we have 72—2, -2 = +2 = — R- 
At the same time, the number of slow electrons increases 
because both the impact electron and the ionization elec- 
tron emerge as slow particles with speed v±\. Because of 
the Coulomb repulsion, the two slow electrons are more 
likely to emerge in opposite directions and we choose the 
rates such that one electron emerges with speed of u_i 
and the other with v+\. 

This choice of model parameters ensures that the en- 
ergy balance during the ionization reaction is not vio- 
lated. As discussed in section III Al the energy of the 
incoming electron is larger that the sum of the energies 
of the escaping electrons because some of the impact en- 
ergy covers the binding energy of the bound electron. In 
the above model, a single ionization reaction transforms 
one electron with speed u+2 into two electrons with, re- 
spectively, speed v+i and speed V-\. The energy of the 
impact electron is 77w+ 2 /2 = 2mAx / At, while the sum 
of the escaping electrons is merely 2mv\/2 = mAx/ At, 
where m is the mass of the electron. So a portion 
mAx/ At of the impact energy covers the binding energy. 
For our model, this is half of the initial impact energy; 
for more general problems other values are possible. 



External Force 



We now derive a discretization of the E^Ji term that 

ov 

models the external force in the Boltzmann equation. We 
start by expanding 



8i_ 

dv 



a v + oxvi 



^m-l^m-l j 



The Galcrkin condition leads to the linear system 



/ m a /3 \ / a \ 

' a P [ ai 

a (3 7 02 

/3 7 a 3 

V /? 7 <5 / \aj 



'J ft!* 

Vvi-Iy 



(12) 



where a = J2 ie s v i>P = E ie5 v t > 7 = E; e s v f and 6 = 



es u i 



vy. 



To calculate the right-hand side of ([12")) we make a 
detour around the continuous representation. We note 
that 



Or 



-t-oc 



l df{x,v,t) 

v z dv ' 

ov 



(13) 



where I £ {0, 1, . . . , m — 1}. Because of our particular 
choice of basis vectors and the fact that there are no 
particles with infinite velocities, we have that 



idf(x,v,t) 

v dv 

Ov 



dv 

f{x,v,t)—dv 



- 3C 

v l f(x,v,t)\±™ 



(14) 



or, in other words, 



df 



dv 



+00 



f(x,v,t)v l 1 dv 



(15) 



With the help of 
( 



iV 





-1 





-1 -1 -1 

— 2w_2 —2v-\ —2vq —2v\ —2vi 

-'W 

3 
1 



-3v_ 2 2 -3w_i 2 -3v 2 Svf 
V -4w_ 2 3 -4w_i 3 -4v 3 -4w 3 



(16) 



-4^ 3 / 



where V = {vo, vi, . . . , v m _i} forms a linear indepen- 
dent set of vectors in M. m . We find the coefficients 
do, ai, . . . , a m -i by enforcing the Galerkin condition. 



V 




we can now define 



V= 



(\ V-2 v- 2 z w-2 3 w-2 4 \ /m a p\ 

'OaO/30 1 
a (3 7 
(3 7 
V P 7 6 J 



1 V„l V-\" V-l" W_l 
2 

2 

\l v 2 v 2 2 



1 Vq V 2 V 3 U 4 
1 Vi Vi 2 Vi 3 V\ A 

V2 3 V2 4 



J 



N, 



In the current paper, we choose a particular set of vec- anc j calculate the external force term as 



tors in V, namely the polynomials {l,v, v 



,m — 1 



} 



discretized in the points w,-. For the five speed example 
the vectors arc, besides their powers of Ax/ At, 



■ *) *! ^ = E(x, t) Vijfrix, t), 



dv 



v = 



1 
1 



(zl\ 



\ 2 / 




1 

UJ 



V 8 J 



f\ 6 \ 


Vl6/ 



where the elements of S denote matrix elements. 



(17) 
(18) 



From cq. (JT5]) , it is clear that we can include the exter- 
nal force as an additional collision term in the right-hand 
side of the lattice Boltzmann equation. 
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D. Flux 

The evolution of the electrical field E(x,t) is deter- 
mined by the net flux j(x,t) of electrons as expressed in 
dU) . We discretize $2$ on a staggered grid with grid points 
halfway between the grid points of the lattice Boltzmann 
model. The flux is defined as the number of particles 
that move between grid points (pass through an inter- 
face) within a single time step. For the five speed model, 
we have 

j(x + Ax/2, t) = f x {x + Ax,t)- f- X {x,t) 

+f 2 (x + Ax,t)-f. 2 (x,t) (19) 
+f 2 (x + 2 Ax, t) - f- 2 {x - Ax, t) 



E. Coupled equations 

The coupled equations |T]) and © for the evolution of 
the electron distribution functions and the electrical field 
is now, after discretization, 

f l {x + v l At,t + At) - fi(x,t) = 

- -(fi(x,t) - f! q (x,t)) +J2^tTZ tJ f J (x,t) 



(E(x - Ax/2, t) + E(x + Ax/2, t)) 



Y^Ai\,,fA.r.t) 



E(x+^,t+At) = E(x+^,t)-Atj(x+^,t), (21) 

where j(x + Ax/2, t) is calculated from (fT9|) . The equa- 
tions are coupled because the electrical field appears as an 
external force in the first equation, while the flux drives 
the evolution of the electrical field in the second equa- 
tion. This coupling makes the evolution of the system 
non-linear. 



IV. A PDE MODEL THROUGH 
CHAPMAN-ENSKOG EXPANSION 

The model (EOJ)— dHU) evolves the electrical field E(x, t) 
and the distribution functions fi£s( x ,t) from t to t + 
At simultaneously. Alternatively, the evolution of the 
distribution functions can be rewritten in terms of the 
corresponding (dimcnsionlcss) velocity moments defined 
as 



Qi{x,t) 



c l ifi(x,t), 



(22) 



where I G {0,1, ... , m— 1}. The zeroth moment Qi=o(x, t) 
corresponds to the electron density p(x, t), i.e. the macro- 



scopic variable of interest. The transformation between 
distribution functions fi and moments Qi can be written 
as a matrix transformation M . In the five speed model, 
this matrix is 



M 



1 1 1 1 1 \ 

-2-10 1 2 1 

4 10 14 
-8-10 1 8 
16 1 1 16/ 



(23) 



such that qi = J2 ie s M nfo and h = TaLq ( M 

An evolution law for gi(x,t) is now easily constructed 
by the following sequence: first transform gi into fi(x,t) 
using M _1 , then use the lattice Boltzmann equation (|2"01) 
to evolve fi(x, t) to fiix, t+At) and, then, transform back 
to the moments Qi(x, t + At). 

It has been observed phenomenological that the ion- 
ization wave can approximately be described by a PDE 
in the density. This suggests that, in practice, the evo- 
lution of these moments is rapidly attracted to a low 
dimensional manifold described by the lowest moment 
Qo(x,t), which is the density. The higher order moments 
have then become functional of this density and the dy- 
namics of the system can effectively be described by the 
evolution of this macroscopic moment. 

In general, however, it is very hard to find analytic ex- 
pressions for this low dimensional description in the form 
of a PDE without making crude approximations. For the 
problem at hand, we illustrate these difficulties in this 
section where we apply the Chapman-Enskog expansion 
and derive a macroscopic PDE in terms of electron den- 
sity. Only after dropping several coupling terms, a closed 
PDE is derived. 

The model discussed in the previous sections can be 
summarized by the lattice Boltzmann equation 

fi(x + CiAx,t + At) -fi(x,t) = 

-- (Mx,t) - f^(x,t)) +Y / A l3 f 3 (x,t) , (24) 

for Vi € S. Here, Aij can be the reaction term i?^, a 
force term Vij , or a combination of both. 

A second order Taylor expansion of the term fo(x + 
CiAx, t + At) in (|24|) around fi(x, t) leads to 

cAx^ + At^ ■ 



+c< Ax At 



dx dt 2 dx 2 

d 2 f % , At 2 d 2 f t 



dxdt 2 dt 2 
= --(/i-/D + E^i' Vie 5. (25) 



jes 



We then expand fi in terms of increasingly higher order 
contributions as follows 



(26) 
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with e a small tracer parameter. In fluid dynamics, e 
typically refers to the Knudsen number. The spatial and 
time derivatives are scaled respectively as 



d_ 

at 



d_ 



d_ 



8U 



and 



d 



d 



dx dxi ' 



(27) 



where we explicitly presume that a zeroth order time 
scale to is present in the system. As we will show later 
on, this scale corresponds to the observed exponential 
growth of the electron density. 

Because of the multiple time scales £o, t\ and t%, all the 
terms in the expansion will couple to all the time scales, 
which complicates the derivation of an effective equation. 

In our search for a reduced second order PDE model, 
we only keep the terms up to second order in e 2 . For the 
same reason, we also drop the second derivative w.r.t. 
time from ([2%]). Substitution of (gHJ) and ([27]) into 
leads to 



eCiAx 
+At 



Of, 



(0) 



-eAt 



dx 

dfj 0) 
dt 



e 2 c;Ax 



Of, 



(i) 



2 c 2 A* 2 d 2 /f 



dx 



dx 2 



91, 



(i) 



dt 



E 2 At 



Of, 



(2) 



dti 



e 2 A£ 



dti 



e'At 



0t 



+eciAxAt 



+e z CiAxAt 

-!(/•»> + 

+ 5>j(/ 

jes 



d 2 f\ 0) 
dxdto 

Q 2f (0) 



iAxAt 



dt 2 
dxdto 



dxdti 



(28) 



We will now group the terms order by order and derive 
expressions for / 4 , fjy and /| and the corresponding 
evolution equations for p at the different time scales. 

We will use the fact that if At and Ax 2 are of the same 
order of magnitude — which is the case for our examples 
— the terms that have factors as At Ax are effectively of 
order Ax 3 and can be neglected compared to terms with 
At or Ax 2 . 



Zeroth order contribution 



and are found as the the solution of the linear system 

£(l-r^)/f =/f. (30) 
jes 

Since f^ q = Wip ijHJ) only depend on the density, we find 
/v = wf p with the weights wf 1 defined by 



t»f ) = E(( 1 -^)" 1 

jes 



(31) 



Since the matrix A can include the ionization reac- 
tion that docs not conserve the particle number, the 
sum of the weights Eies w i ^ s n0 ^ necessarily equal 
to one. This means that J2ies H ^ P- We choose 
to rescale the weights u)( with a normalization factor 



= E W65 ((i-tA) 



such that 



J2ies fi = P- With this rescaling we find the zeroth 
order term of the Chapman-Enskog expansion 

/r=«f ) P=^E(( 1 -^)" 1 ). W TP (32) 



jes 



The rescaling, however, forces us to reconsider equa- 
tion ([30]) because , as defined above, fails to be a 
solution. Still, we keep our f- ' of ([52")) as our zeroth- 
order term and find for the evolution of the system at 
time scale to 



At 



Of, 



Co) 



dt 



jes 



jes fees ' 

where the growth factor is 



(33) 



a= (l-l/A0/(rAt). 



(34) 



From expansion ([28]) , we collect the zeroth order terms 

1 ^(0) ,e<A , A AO) 



At- 



dt 



(29) 



We choose ff 1 such that the right-hand side of ([2"4"| is 
zero; these distributions will not evolve on time scale to 



Summation of (|3"3")l over the set S leads to the zeroth 
order PDE for the evolution of p 



dp_ 

dt 



ap. 



(35) 



This is a growth equation if a is positive, which is the 
case for the ionization reaction. 
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2. First order contribution 



3. Second order contribution 



To derive the first order equation, we collect the terms 
that arc first order in e from (1281) 



;A. 



Co) 



dx 



At 



dfl 



(i) 



,(o) 



df> 

-At-^r h a At Ax 



dt 

dxdt 
(i) 



(36) 



We drop the term CiAtAxd 2 ff^/{dxdto) because it is of 
order AtAx, which is smaller than the other terms that 
are first order in At or Ax. The second term we neglect 
is Atdfi /dto because we will show below, a postiori, 
that it is also of order AtAx. 



We now have 

CiAx^r— + At ■ 



dx dti \ t 

jes 



F (l) 



what leads to a first order term 



/<"=E(MA^r') (j ( 

jes J \ 



<9/j 0) dL 
Cj Ax^- + At^f 
dx dti 



(o)' 



(37) 

We now see that it is justified to neglect the term 
Atdf^' /dto in ([36]) because it is of order AtAx. Using 
([32]) and J2i^s fi^ ~ *-* ( tnc l attcr because X^gs /j = 



E, e5 (/f +e/f ) +e 2 /fj = P ^d £ ieS /r = p), 
we find the following PDE at time scale t\ for the evolu- 
tion of the system 



p(0) 



dti dx 
where the advection coefficient c equals 



(38) 



C = 



^M^I^X^Ax 
a^((-l/r + ^)-%«,r Af 



Finally, we derive the second order evolution. We col- 
lect from (1281 the second order terms and find 



<9f- (1) c 2 Ax 2 d 2 f (0) df {2) 
Ci Ax^4- + 1 A + At h 



dx 2 dx 2 



dto 



df (D d AO) d 2AD 

+At^— + At-£- + Cl AtA Xl -±- 
dti dt2 dxdto 

r) 2 f (0) 1 
^AtAx^ = --£(l-^,)/f (41) 



3 £5 



The terms d 2 / (dxdt ), d 2 / {dxdh) , and 

df- 1] /dh (when replacing /f ) by jST}) are of order 
AtAx, which is smaller than Ax 2 for our parameter 

(2) 

settings. We also neglect Atdf\ /dto because it can be 
shown, again a postiori, that it is of order AtAx. 

The second order expansion term then becomes 
df {1) c 2 Ax 2 d 2 f {0) 



dx ' 2 <9x 2 
aA°) 1 



When replacing /j and fjy by (|40f and (|32|) . we get 



c 4 Ax £ ((-1/r + A)" 1 ) t (c,- Ax - CAt) 

i c|Ax 2 Wl (0) d 2 p Ai (0)9p_ 
+ 2 dx 2 + 1 dt 2 



= -^(l-r^)/f(43) 



and the expression for the second order term becomes 



fP = E s ^ E ^ fe ( cfcAa; - CAt ) A 



j£5 jeS 



(o) d 2 p 
dx 2 



— ' ^ 2 dx 2 



(44) 



where we use Bjj = ((— f/r + A) If we define 



With the help of (|38[) . the first order contribution (|371 
can be written alternatively as 

if = E ((-Vt + A)" 1 ) «,'°>(c, Ax - CAt)|. 



J£5 



(40) 



p= E B ijCj B jk (c k Ax-CAt)A 



xw 



(0) 



(45) 



J2 B^vf Ax 2 /2 / - E ^< )A * I > 
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we obtain the evolution at time scale ti (because 

f(2) 



0.05 



0) 



dp 



V 



d 2 x ' 



(46) 



where T> acts as a diffusion coefficient. 

We are now in the position to combine the evolution at 
the different timcscales to, t\ and £2 into a single PDE. 
Because the matrix A of the model equation (|24|) con- 
tains both the reaction term and the external force, the 
transport coefficients a, C and T> will depend on the local 
electrical field E(x,t). For our example the dependence 
of the transport coefficients is shown in figure^ we find 
that V hardly depends on the strength of the field and 
can be set equal to the electron diffusion D used in © 
to define the relaxation of the lattice Boltzmann model. 
In a similar way, we find that c, the transport coefficient 
of the advection term, is approximately equal to the —E, 
the local electrical field that causes the drift. Only the 
growth factor a, defined in (|34p . depends on the strength 
of the local electrical field. With the help of these obser- 
vations we get the coupled PDE 



dp 

Of 
dE 

~dt 



a(E(x,t))p + E(x,t) 



dp ^d 2 p 



O.r 



D- 



dx 2 



-E(x,t)p- D 



dp 
Of 



(47) 



for the evolution of E{x, t) and p{x, t). The second equa- 
tion expresses the flux with the help of the transport 
coefficients. 

These coupled equations arc similar to minimal 
streamer model of Ebert, Van Saarloos and Caroli 
except that the growth rate is now defined by (|34|) in- 
stead of the Townsend approximation. In figure [H we 
illustrate how the growth coefficient depends on the lo- 
cal electrical field and compare with a Townsend ap- 
proximation. We find that a Townsend reaction term 
0.111 ■ \E\ exp(—l/\E\) approximately describes a similar 
growth term as the PDE model derived from the lattice 
Boltzmann model. 



4- Traveling wave solutions 
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Figure 2: Top: The growth factor a(E) (solid) of the 
PDE model derived from the lattice Boltzmann model. The 
growth depends on the strength of the local electrical field 
and is similar to the Townsend approximation with 0.111 ■ 
\E\ exp(— l/|_E|)(dashed). We have a reaction rate R = 100 
and model parameters given in section IVHI Middle: The ad- 
vection coefficient C is equal to the external field — E. Bottom: 
The diffusion coefficient T> only changes with the external field 
in the fourth significant figure. 



The system (|47[) is non-linear and it is well-known 
that it has a one-parameter family of front solutions that 
translate uniformly with a speed c 0] • There is a minimal 
speed c* that is usually found by looking at the asymp- 
totic region — > +00. In this limit, the electrical field 
becomes constant and is denoted by E + — the same 
notation as in [l[ — and the equation for the electron 
density becomes 



where the transport coefficients have become constants. 
In a co-moving coordinate frame that travels with the 
same speed c along the x-axis, we define £ = x — ct. 
The PDE (|4"5)) becomes stationary and the solution in 
the asymptotic region fits 

= a(^)p(0 + (c + ^)^ + U^. (49) 
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The latter is a second order ODE that can be transformed 
into two coupled first order ODEs by denoting dp/d£ as 
v and p as u . The system of coupled equations is 



Algorithm 1 Constrained runs scheme for LBM 







1 

C+E-* 
D 



(50) 



where and denote derivatives of, respectively, u and 
V to £. The matrix has two eigenvalues 



A± 



-(c + E+) ± yj{c + E+Y -ADa{E+) 
2D ' 



There are two cases, if (c + E + ) 2 < ADa(E + ) the two 
eigenvalues are complex, otherwise, they are real. 

The electron density in the asymptotic region is now a 
linear combination of two exponentials 



lim p(x) = Ae x + X + Be > 

x— »+oo 



(51) 



When the two eigenvalues are complex, the asymptotic 
density is oscillating and can becomes negative. This is 
unphysical because we cannot have a negative number of 
particles and it is concluded that the speed c has to be 
above a minimal speed 



> c* = c(E+) + 2y/Da{E+), 



(52) 



to keep both eigenvalues real. Note that both eigenvalues 
A+ and A_ coalesce at the critical speed c = c* . 



V. THE COARSE-GRAINED TIME-STEPPER 



In this section, we describe an alternative way to per- 
form the analysis of the macroscopic behavior of the sys- 
tem. It is based on the work of Kevrekidis et al. Q who 
developed a coarse-grained time stepper (CGTS), which 
is an effective evolution law for the density. This evolu- 
tion law J- is not an analytic expression such as a PDE, 
but the following sequence of computational steps: 1) 
lifting, 2) simulation and 3) restriction, denoted by the 
operators p, LBM and M. respectively (figure [3]). Note 
that the simulation time AT is in general a multiple of 
At, the lattice Boltzmann time step. Formally, this is 
written as 



U(x,t + AT) = T(U(x, t), AT) 

= M(LBM(p(U(x,t)),AT)), 



(53) 



where we have introduced U(x,t) = {p(x,t),E{x,t)) as 
a shorthand notation. The time-stepper T evolves the 
macroscopic density p(x,t) = X^e<s fi( x ->t) an d the elec- 
trical field E(x, t) from time t to t + AT. 



initialize /j ' = w^ q p(x, t) 
repeat 



Vi e S 



flk+i] = LBM(/W) a single LBM step 
gl k +i-} — Mf[b+i\ map into moments 



p{x,t) 



reset the density 



flk+M — m 1 g\- k+1 ' map into distributions 
until convergence heuristic 

Table I: Lifting. The constrained runs algorithm computes 
the distribution functions fi(x,t) corresponding to a given 
density p(x, i). The superscript k indicates the iteration num- 
ber. 



A. Lifting 

Since the electrical field E(x, t) is the same in both the 
lattice Boltzmann and the macroscopic model, we can 
ignore it for the discussion of the lifting and restriction 
operators. In the lifting step, the particle distribution 
functions are initialized starting from the initial density 



/' 



p(Xj,t) H-> fi(Xj,t) 



with i £ S, m the number of speeds in S, and j £ 
{1,2,..., n} denoting the discrete spatial grid points. Be- 
cause lifting is a one-to-many mapping problem, it is the 
most critical step in the coarse-grained time-stepper. We 
use the constrained runs scheme, an algorithm proposed 
in in the context of singularly perturbed systems. 
Here, it is wrapped around a single time step At of the 
lattice Boltzmann model [4|. 

The procedure is given in table [H Starting from an ini- 
tial guess p(x, t) for the density, we obtain initial guesses 
for the distribution functions using the BGK equilibrium 
weights ([HI). This choice determines the initial guess for 
the higher order moments through the transformation 
matrix M, see equation (|2"3"|) . (In principle, the initial 
guesses for the higher order moments can be chosen ar- 
bitrarily; the scheme is designed precisely to converge to 
the correct value of these moments for the given density. ) 

We then perform the following iteration. First, we use 
the lattice Boltzmann model to evolve f$ s from t to 
t + At. The result is transformed back into the moment 
representation by a matrix multiplication with M , which 
gives us g[ fc+1 l. Next, the zeroth moment of the vector 
i s reset to the initial value p(x,t). Transforming 
this modified moment vector gl fc+1 l back into distribu- 
tion functions gives us the next /jgj • We repeat this 
iteration until the higher order moments have converged. 

The convergence behavior of the constrained runs algo- 
rithm from table U is analyzed in Q for one-dimensional 
reaction-diffusion lattice Boltzmann models with S = 
{ — 1,0,1} (D1Q3 stencil) and a density dependent re- 
action term. For such systems, the algorithm is uncondi- 
tionally stable and converges up to the first order terms in 
the Chapman-Enskog expansion of the distribution func- 
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Figure 3: A summary of the different steps involved in finding the traveling wave solutions as fixed points of the coarse-grained 
time-stepper. We start with an initial guess for the density in the left bottom (panel a). This density is mapped to the 
corresponding components of the one-particle distribution function using the constrained runs lifting operator. The initial 
conditions (panel b) are then evolved with the full lattice Boltzmann model over a time AT. Each component travels over a 
distance cAT and arrives at panel c. In the next step, the density at time t + AT is extracted using the restriction operator 
(panel d). The resulting density is shifted back over a distance cAT to arrive at the original position. The traveling wave 
solution should be invariant under this sequence of operations and is formulated as a fixed point. 



tions. The convergence rate is 1 — l/r|, i.e. the same 
rate at which the lattice Boltzmann simulation relaxes 
towards the diffusive BGK equilibrium. 

Below, we extend the results from Q to the current 
five speed model with the velocity dependent reaction 
term (fTTj) . In the absence of an electrical force field, the 
distributions f±2 for the fast particles evolve as (|2T))) . i.e. 

f ±2 (x, ± 2A.t, t + At) = (1 AtR)f ±2 ( Xj ,t) 

T 

+ -w^ 2 p(xj,t), (54) 
r 

where the second term is "frozen" because in each iter- 
ation of the constrained runs algorithm, the density is 
reset to its original value. Because the LBM propaga- 
tion of distributions is a conservative operation [16| , this 
iteration is linearly stable if 

II AtR\ < 1. (55) 

r 

The distributions f±\ evolve as 

f ±1 ( Xj ± Ax,t + At) = (l--)f ±1 ( Xj ,t) 

T 

+ -w% lP ( Xj ,t) + AtR(f +2 { Xj ,t) +/_ 2 (ar i ,t)),(56) 

where the number of slow particles is increased pro- 
portional to the number of fast particles because of 



the ionization reaction. Again, the density value in 
the BGK equilibrium in (|5^|) is "frozen" to the initial 
value. Because the convergence rate for the f±2 compo- 
nents is given by (|55p . equation (|56[) converges at a rate 
|1 — 1/r — AtR\ if this value dominates over |1 — l/r|, or 
at a rate |1 — l/r| if the latter is dominant over (|55|) . 

So far, we have no formal proof for the convergence 
when an electrical field is present, but we can illustrate 
the convergence of the algorithm for the full system (j2"0)) 
(f2"Tj) numerically for the parameter settings from sec- 
tion [Vni The figure is produced as follows. We first ex- 
tract the velocity moments (|22[) from a full lattice Boltz- 
mann simulation that has evolved from an initial state 
for several thousand time steps. Subsequently, we use 
the obtained density p as the initial condition for an- 
other lattice Boltzmann simulation and use algorithm U 
for its initialization; the distribution functions of the orig- 
inal simulation are considered to be the "exact" solution. 
Figure [4] (top) plots the norm of the error between the 
constrained runs state and the state of the first lattice 
Boltzmann simulation that it tries to reconstruct. We 
also observe that after initialization, when evolving the 
full lattice Boltzmann system both from the "exact" and 
the re-initialized distributions, the error between the first 
and second system further decreases (figure 21 bottom). 
The same observation was made in Q . Note that we have 
no analytical expression for the initial state returned by 
the constrained runs scheme in this setting, but we be- 
lieve that the results on the accuracy from [J] generalize 
and that the obtained initial state is a first order approx- 
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Figure 4: Convergence of the constrained runs algorithm for 
the five speed ionization model both during lifting and simula- 
tion step. Top: Error (2-norm) in the lifted distribution func- 
tions (circles) and the flux (squares) as a function of the num- 
ber constrained runs iterations. After an initial convergence 
with rate |1 — 1/t — AtR\, the error stagnates after approx- 
imately 25 iterations. Bottom: Difference (2-norm) between 
the distribution functions (circles) and the flux (squares) of 
the lattice Boltzmann simulation that started from the stag- 
nation point of the top figure and the original simulation. 
Again the error decreases for approximately 25 lattice Boltz- 
mann steps. Note however that after stagnation there is a 
slow evolution because the macroscopic fields also evolve. 



Restriction 



In the last step, the restriction step, we extract the 
macroscopic variables from the result of the simulation. 
The macroscopic density at time t + AT is then 



M 



-> R n : ft (xj , t + AT) h-> 
p(x j ,t + AT) = J2fi(x j ,t + AT). (57) 



ies 



VI. 



THE TRAVELING WAVES AS A FIXED 
POINT PROBLEM 



In this section we describe the methodology outlined 
in Q to find the traveling wave solutions of the coarse- 
grained time-stepper f(U(x,t)) defined in sectionfVj If 
a traveling wave solution with a speed c of !F(U(x, t)) is 
evolved over a time AT, the solution has shifted over a 
distance cAT. We define a shift-back operator that 
shifts the solution back over a distance ijj. 

o-^ : U(x,t) h-> a^(U(x,t)) = U(x,t) + ipd x U(x,t), 

where we implement the shift-back by using the charac- 
teristic solution of dtU(x, t) + tpd x U(x, t) in the forward 
Euler time discretization. 

This shift-back operator is combined with the coarse- 
grained time-stepper to write a non-linear system for the 
traveling wave in the co-moving coordinate system with 
£ = x — ct 



U(0-<J cAT CF(tf(0,AT)) = 0. 



(58) 



imation of the Chapman-Enskog relations. 



B. Simulation 

In the simulation step, the initial distributions ob- 
tained from the lifting step, are evolved for a coarse- 
grained macroscopic time step AT using the lattice 
Boltzmann model discussed in the previous sections. 
This step is denoted as 

fi( Xj ,t + AT) = LBM(fi{ Xj , t), AT), 

where the number of time steps depends on the ratio 
of AT/ At. When the ratio is not an integer, a linear 
interpolation is used between subsequent steps. 



This equation expresses the sequence of computational 
steps as illustrated in figure [3] This system, however, is 
singular because any translate of a solution will also be 
a solution of (J3HJ> @j- 

To get a regular system we add phase (pinning) con- 
dition p(U) and a rcgularization parameter a as an ad- 
ditional unknown, as discussed in [3| . The resulting non- 
linear system is 

where the phase condition p(U) is defined as 

p(u) = / u(Z)dtU n f(Z)dZ. 

This condition minimizes phase shifts with respect to the 
reference solution /7 re /(£). 
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A. Preconditioned Newton-GMRES 

We solve the non-linear system (|59[) using a Newton- 
ftaphson method, 

f = i/M + dlfW 



where the corrections dU^ and da^ are calculated by 
solving, in each Newton iteration, a linear system of the 
form 



I- J(C/W,a[ fe l) d 6 U^\ 
-G(U [k \a [k] ). 



dUW \ 
daW ) 



(60) 



This system is the linearization of G(U, a) around the 
point ({/[ fe l, a^) and J{U^ k \a^) denotes the Jacobian 
of (t c &t(F(U^ , a'")). Since T is defined as a sequence 
of computational steps, it is impossible to construct the 
Jacobian J explicitly. We therefore use a Krylov method 
(GMRES) that only requires its application to a vector 
v, which can be estimated as 



(I-J(U,AT))vn 

<JcAT (HU + ev, AT)) - a cAT {T{U, AT)) 



(61) 



Since the convergence rate of GMRES depends sen- 
sitively on the spectral properties of the system matrix 
(|6U)) . we propose to precondition the linear system ([ST)]) 
with a rough macroscopic model based on a PDE to speed 
up the convergence [3j. In section HVi we derived an 
approximate PDE model using a Chapman-Enskog ex- 
pansion. We define a time-stepper for this approximate 
model as F(U(x, i), AT), we can again write a non-linear 
system 



G(U, a) 



U-a cAT {F(U,AT)) = 
P (U) = 



(62) 



in which we have replaced the coarse-grained time- 
stepper by a time-stepper for the approximate PDE 
model. The solution of this system will look very similar 
to the solution of the full model, but will differ in places 
where the approximations made during the derivation of 
the PDE model fail. The linearization of this problem 
leads to a matrix problem for the Newton corrections 
dU [k] and da [k] as in (|60[) with very similar spectral prop- 
erties. The Jacobian however, is now known analytically 
and can be inverted easily. The idea is to use this matrix 
as a preconditioner of the linear system that is solved 
each Newton iteration. The preconditioned system reads 

(^W,aW))' 1 A(^W | aW)(^J) 

= - (7\/(C/[ fc ],a[ fc ])) _1 G(J7 [fcl ,a [fcl ), (63) 
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Figure 5: The traveling wave solution for c — 1.3 and R = 100. 
It is a fixed point of sequential application of the evolution 
with the coarse-grained time-stepper over a time AT and the 
shift-back over a distance cAT. 



where (M(U^ k \a^)) is the inverse of the matrix of 
the linearization of (f6"2"| and A^^, w k ') denotes the 
linear system of ([6H|) . Because the spectral proper- 
ties of A(Z/t*],aW) and M{U^ k \a^) are so similar that 

(M(i/[ fc ],a[ fc ])) -1 A(UW,a.W) has spectral properties fa- 
vorable for the convergence of GMRES. Detailed numer- 
ical experiments, showing the spectral properties of the 
linear systems and the GMRES convergence, are reported 
in @. 
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B. The minimal speed and the coarse-grained 
time-stepper 

In section IV, we derived the minimal speed c* (|52|) 
of the uniformly translating front solution of the PDE 
model in terms of the asymptotic transport coefficients. 
Solutions with a speed below this critical value have a 
density that oscillates in the asymptotic region, which 
leads to unphysical negative densities. 

For the coarse-grained time-stepper, the asymptotic 
transport coefficients are not available and no analytic 
expression for the minimal speed can be found. We pro- 
pose to use the time-stepper and its fixed point solutions 
to determine the minimal speed. We vary the imposed 
speed of the shift-back operator <7 c at and monitor the 
solutions of the fixed point problem in the asymptotic 
region. As the imposed speed c falls below the minimal 
speed c* the solutions will become oscillatory because 
the solution in the asymptotic region is a combination of 
two exponentials (|5ip . whose exponents X± coalesce and 
become complex at the critical speed. 

We can extract the two exponents from the solution in 
the asymptotic region, if we assume that it fits a second 
order ODE of the form 

d 2 p dp 
—— = aip + a 2 — . 

ox z ox 

The coefficients a\ and a 2 are found by taking the fixed 
point solution in two grid points, where we estimate the 
spatial derivatives using finite differences. This allows us 
to formulate a 2-by-2 system for a\ and a 2 . 
The eigenvalues of the 2-by-2 matrix 

( ° 1 ) (64) 

arc then A + and A_, which coalesce at the critical speed. 

For a speed c far above the critical speed c* , this 
method estimates only one of the two exponents with 
confidence. Indeed, above the critical speed the solution 
in the asymptotic region is a combination of two decay- 
ing exponentials with different exponents. However, one 
of them is slowly decaying, while the other decays fast. 
Far away from the critical speed, the method will only 
detect the slowly decaying exponential and a fit with a 
first order ODE would be sufficient. Near the critical 
speed, however, the solution is a combination of both ex- 
ponentials that decay with comparable rates and both 
eigenvalues can be reliably extracted from the solution. 

VII. NUMERICAL RESULTS 

As an illustration, we look at a one-dimensional lattice 
Boltzmann model on a grid with TV = 1600, grid dis- 
tance Ax = 0.4 and a time step At = 0.008. We look at 
a model with five velocities with (S = {—2, — 1, 0, — 1, 2}), 
weights w\ q = {0, 1/4, 2/4, 1/4, 0} and an electron diffu- 



sion coefficient of D = 1.0. This leads to a relaxation 
parameter of r = 0.8 or u> = 1.25. Note that with this 
choice of equilibrium weights only slow particles exist in 
the absence of external fields. 

We enforce boundary conditions at the level of the lat- 
tice Boltzmann model where we use homogeneous Dirich- 
let boundaries at the right and no-flux boundaries at the 
left. The electrical field is kept constant at E + = — 1.0 
at the right and at the left we require that d 2 E/dx 2 = 0. 

We first compute the traveling wave with speed c = 
1.30 for a reaction rate R = 100, which is shown in figure 
[5] As an initial guess for the Newton procedure, we take 

p{x,t) = 0.025/(1 + exp(0.15(x-L 2/3))) , , 
E{x,t) = -1/(1 + exp(0.05(x-L 5/9))) [ ' 

To assess the overall performance of the method, we com- 
pute the total number of required lattice Boltzmann time 
steps. We observe that about 5 Newton steps are needed. 
Each Newton step, in turn, requires the solution of a 
linear system with the help of a preconditioned Krylov 
subspacc. On average about 40 GMRES iterations are 
required to solve the linear system with a tolerance of 
1 • 10~ 12 . Each GMRES iteration requires an evaluation 
of the CGTS, which costs about 50 lattice Boltzmann 
iterations — 25 for the lifting and 25 for the simulation. 
This leads to a total of 10 000 evolutions with the lattice 
Boltzmann system to find a single fixed point. Detailed 
figures about convergence are given in Q. 

Next, we vary the cross section R, which is related to 
the number of times an ionization reaction appears, and 
study the effect on the critical velocity of the traveling 
wave solution. We increase R from 60 to 100, which cor- 
responds at the level of the PDE to an increase of the 
growth transport coefficient from a(E + ) = 0.02015 to 
a(E+) = 0.03707, with E+ = -1.0. For this range of re- 
action rates with the chosen value of At, the constrained 
runs algorithm always converges because all the eigenval- 
ues of the Jacobian are smaller than 1 1 — 1/r — AtR\ < 1 , 
as discussed in section IV Al 

We now turn to the numerical computation of the criti- 
cal velocity with the method proposed in section lVI Bl In 
tabic [Til we show the critical speed determined for a se- 
ries of ionization strengths between R = 60 and R = 100. 
As outlined in section lVI Bl the critical velocity for each 
value of R is found by performing a numerical contin- 
uation with the wave speed c as a free parameter and 
monitoring the eigenvalues of the 2-by-2 matrix (|64p . 
For comparison, we also show the minimal speed as ob- 
tained using the approximate analytical expressions for 
the transport coefficients, which we found through the 
Chapman-Enskog expansion. We observe a small differ- 
ence between the two results. The critical speed that 
results from the Chapman-Enskog expansion is higher 
than the critical speed that is computed using the CGTS. 
Note that both methods make approximations and it is 
not clear which method is more exact. 

A short discussion about the accuracy of the shift- 
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R 


c* with 


c* with 




Chapman-Enskog 


CGTS 


60 


1.3351 


1.3177 


70 


1.3496 


1.3318 


80 


1.3609 


1.3433 


90 


1.3696 


1.3517 


100 


1.3755 


1.3566 



Table II: The critical speed for different values of the ioniza- 
tion rate R. 

back operator is necessary. Since we have implemented 
the shift-back operator using a forward Euler approx- 
imation of the characteristic solutions of the equation 
ut + cu x = 0, the error in the critical speed grows linearly 
with the time step AT. Therefore, we should take AT 
as small as possible, i.e. the minimal possible number 
of lattice Boltzmann steps in the simulation step. (Note 
that accuracy and efficiency go together here.) However, 
we cannot take less than 25 LBM steps because these 
steps are necessary to reduce the error in the lifting (see 
figure . 

VIII. DISCUSSION AND CONCLUSIONS 

In one dimension, an initial seed of electrons in a strong 
electrical field will evolve into a streamer front that trav- 
els with a constant speed c. Before the front, the density 
is zero and the electrical field is constant; behind the 
front the electrical field is shielded and there is a sur- 
plus of electrons. In this article we extended the min- 
imal streamer model of Ebert et al. [l[ to add some 
details about the microscopic physics. To this end, we 
replaced the reaction-diffusion PDE with a lattice Boltz- 
mann model that has a velocity dependent reaction term. 
The reaction rates are a chosen to model the ionization 
reaction, where fast particles have a given probability to 
undergo an ionization reaction and create two slow par- 
ticles. The electrical field changes simultaneously as a 
result of the charge creation. 

This macroscopic behavior of the model was analyzed 
with two methods. The first is the more traditional anal- 
ysis based on a Chapman-Enskog expansion that derives 
an approximate PDE model and the corresponding trans- 
port coefficients. The resulting PDE model is very simi- 
lar to the minimal streamer model of Ebert, Van Saarloos 
and Caroli, where the ionization rate depends on the lo- 
cal electrical field. Based on the transport coefficients, we 
found the dependence of the critical speed (below which 
the traveling waves are unphysical) and its dependence 
on the strength of the ionization cross section. 

Our second method is a computational method based 
on the coarse-grained time-stepper, which defines the 
effective evolution law as a sequence of computational 
steps. The traveling wave solutions were formulated as 
fixed points of this coarse-grained time-stepper combined 



with a shift-back operator. By varying the applied shift- 
back, we again found the critical speed, and we demon- 
strated how to calculate its dependence on the cross sec- 
tion strength numerically. 

We showed that the coarse-grained time-stepper pro- 
vides a viable alternative to the traditional Chapman- 
Enskog analysis. However, in this paper, we did not 
make any statements about which of the two methods 
provides the most accurate information. Both methods 
make approximations to obtain the critical velocity of 
the proposed lattice Boltzmann model. In contrast to 
the theoretical approximations in the Chapman-Enskog 
analysis, these approximations are due to numerical ac- 
curacy for the coarse-grained time-stepper. 

Further progress in the accuracy can be made in several 
ways. First, the lifting can be done more accurately if a 
higher order constrained runs algorithm is used. Such a 
method would take multiple steps with the lattice Boltz- 
mann model and uses these multiplepoints to provide us 
with an approximate initial state [17] . We expect that 
these higher order lifting procedures will reproduce more 
than two terms in the Chapman-Enskog series and there- 
fore provide a better initial state for the simulation step. 
This will allow limiting the number of simulation steps, 
which will immediately improve the accuracy of the com- 
puted solutions and the critical speed. 

Although the model problem in this paper is non- 
trivial, and the Chapman-Enskog expansion is tedious, 
we emphasize that our method is mainly developed with 
applications in sight where the traditional Chapman- 
Enskog analysis is completely intractable. The model 
problem in this paper both allows us to analyze the pro- 
posed methods and provide directions for further devel- 
opment. In a forthcoming article, we will do a complete 
analysis of the traveling wave solutions in the proposed 
model and illustrate more extensively how the macro- 
scopic behavior depends on the parameters of the micro- 
scopic model. In our calculations, we have also observed 
the positively charged fronts that move in the opposite 
direction of the electrical field. These solutions have been 
discussed by Ebert et al in [l[ and our methods are also 
able to locate these states. 

We expect that most of the results presented in this 
paper will remain valid when the number of velocities in 
the discretization of the Boltzmann model is increased. 
With additional velocities, it is possible to model the 
cross sections of colliding particles with more detail and 
study other reactions than the ionization reaction. It 
would also be interesting to study models with multiple 
species where the collision rates between the particles 
are velocity dependent. This would allow us to include 
photo-ionization effects into the minimal streamer model, 
an important effect that is neglected in the current model 
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